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Abstract 

We present predietions for double-quarkonium produetion in the kinematieal region relevant 
for the proposed fixed-target experiment using the LHC beams (dubbed as AFTER@LHC). These 
inelude all spin-triplet 5-wave eharmonium andbottomoniumpairs, i.e. \J/{n\S) + if/{n 2 S), if/{n\S)-<r 
T(mi5) and T(mi5) -l- T(m25) with ni,n 2 = 1,2 and mi,m 2 = 1,2,3. We ealeulate the eontribu- 
tions from double-parton seatterings and single-parton scatterings. With an integrated luminosity 
of 20 fb“^ to be collected at AFTER@LHC, we find that the yields for double-charmonium pro¬ 
duction are large enough for differential distribution measurements. We discuss some differential 
distributions for J/ij/ -l- J/iJ/ production, which can help to study the physics of double-parton and 
single-parton scatterings in a new energy range and which might also be sensitive to double intrin¬ 
sic cc coalescence at large negative Feynman x. 

Keywords: Quarkonium production, Double-parton scattering, QCD 
PACS: 12.38.Bx, 14.40.Gx, 13.85.Ni 


1. Introduction 

Heavy-quarkonium production is typically a multi-scale process, which involves both short- 
and long-distance facets of the strong interaction. This particularity makes heavy-quarkonium 
production an ideal probe to study Quantum Chromodynamics (QCD) in its perturbative and non- 
perturbative regimes simultaneously. Studies have extensively been performed at collider and 
fixed-target energies in proton-proton, proton-nucleus and nucleus-nucleus collisions (see reviews 
e.g. Refs. [1, 2, 3]). The associated production of heavy quarkonium is a very interesting process 
not only because it provides a way to pin down the heavy-quarkonium production mechanism 
but also because it can help to understand a new dynamics of hadron collisions appearing at high 
energies, where multiple scatterings of partons (MPS) happen simultaneously, among which the 
most likely is of course two short-distance interactions from a single hadron-hadron collision - 
double-parton scattering (DPS). A number of experimental studies relevant for DPS analyses with 
heavy quarkonia have recently been carried out such as 7/i/^ -l- W [4], Jjifj + Z\5], 7/i/^-l-charm [6] 
and Jjifj + Jlif/ \1] production. 
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In particular, the latter process, i.e. double-quarkonium production, is of specific interest. It 
provides an original tool to study the quarkonium production from the conventional single-parton 
scatterings (SPSs), whose contribution has theoretically been studied in many works [8, 9, 10, 11, 
12,13,14,15,16,17,18, 19]. Moreover, it has been claimed in Refs. [20,21,22,23,24,25,18, 19] 
that DPS contributions should be a significant source of Jjif/ + Jjif/, especially at high energies 
where there is a high gluon flux. On the experimental side, the spin-triplet S -waves {e.g. J/if/, if/', 
T(nS)) provide clean signatures with their small background when they are studied in their decay 
into muon pairs. They are easy to trigger on, in contrast to hadronic jets and open-charm meson 
productions, which require either good calorimetry or good particle identification. 

A first comprehensive comparison between experiments [26, 7, 27] and theory for 7/i/r-pair 
production at the Tevatron and the LHC has been performed in Ref. [18], where we have pointed 
out that this observable could be used to probe different mechanisms in different kinematical re¬ 
gions. We noted that the direct DPS measurement by DO collaboration [7] -looking at the rapidity- 
difference spectrum- is consistent with the y/i/^-pair measurement by the CMS collaboration [27] 
and, as we will discuss later on, compatible with rather large DPS rates. On the other hand, as 
we advocated in [16], one cannot draw a definite conclusion on the presence of DPS in the early 
LHCb data [26] with their relatively low statistics. 

In this context, we find it important to study the potentialities offered by the use of the 7 TeV 
proton LHC beams in the fixed-target mode to study quarkonium-pair production. Its multi-TeV 
beams indeed allow one to study p + p, p + d and p + A collisions at a centre-of-mass energy 
■y^SNN - 115 GeV as well as Pb -l- p and Pb -l- A collisions at -^^Snn - 72 GeV, with the high 
precision typical of the fixed-target mode. It has indeed been advocated in [28, 29] that such a 
facility, referred to as AFTER®LHC, would become a quarkonium, prompt photon and heavy- 
flavour observatory thanks to its large expected luminosity (for recent phenomenological studies, 
see [30, 31, 32, 33, 34, 35, 36, 37, 38, 39]). A first feasibility study for quarkonium production 
was presented in [40] and demonstrated that a LHCb-like detector would perform extremely well 
in the fixed-target mode. Similar performances are expected for quarkonium-pair production. 

Integrated luminosities as large as 20 fb“' [28] can be delivered during a one-year run of 
p -I- H collisions with a bent crystal to extract the beam [41]. The LHC beam can also go through 
an internal-gas-target systemk Conservatively sticking to gas pressures already reachable now, 
yearly integrated luminosities reach 100 pb“^ With a designed target cell similar to that of HER¬ 
MES [45], a few fb“^ yr“' are probably also easily reachable [46]. We have reported in Tab. 1 the 
instantaneous and yearly integrated luminosities expected with the proton beams on various target 
species of various thicknesses, for both options. 

The structure of this paper is as follows. In section 2, we detail and justify our methodology 
to compute both DPS and SPS contributions to quarkonium-pair production. Section 3 contains a 
general discussion of the interest to look at DPS vs SPS contributions at different energies. Section 
4 presents a comparison between results up to a] and a^. This prepares the discussion of our results 
at = 115 GeV relevant for APTER@EHC in Section 5. Section 6 gathers our conclusions. 


'This is in fact already tested at low gas pressures by the LHCb collaboration in order to monitor the luminosity 
of the beam [42, 43, 44]. 
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Beam 

Target 

Thickness 

(cm) 

P 

(g.cm-3) 

£ 

(yUb-^S-l) 

(pb-'.y-i) 

P 

Liquid H 

100 

0.068 

2000 

20000 

Beam 

Target 

Usable gas zone 
(cm) 

Pressure 

(Bar) 

£ 

Oub-^s-i) 

JX 

(pb-'.y-‘) 

P 

perfect gas 

100 

io-‘^ 

10 

100 


Table 1; Expected luminosities obtained for a 7 TeV proton beam extracted by means of a bent crystal or obtained 
with an internal gas target with a pressure similar to that of SMOG@LHCb [43]. 


2. Methodology 

In this section, we explain the main ingredients used to compute the rates for double-quarkonium 
production at AFTER@LHC, which closely follows from our previous work in Ref. [18]. 

2.1. Double-parton scatterings 

The description of such a mechanism is usually done by assuming that DPSs can be factorised 
into two single-parton scatterings (SPS) resulting each in the production of a quarkonium. This 
can be seen as a first rough approximation which can however be justified by the fact that possible 
unfactorisable corrections due to parton correlations could be small at small x. In the case of the 
double-quarkonium production, the master formula from which one starts under the factorisation 
assumption is (see e.g. Ref. [24]) 

(Tq^q^ = -I dx\dx2dx\dx2cp'\)\(f'\)2cf'\i 

X Tij{xi,X2, bi, b 2 ) &‘f^{xi,x\)a^l{x2, x'^) T^/Cx], 4,bi - b, b 2 - b), (1) 

where Yij{x\,X 2 , bi, b 2 ) is the generalised double distributions with the longitudinal fractions X\,X 2 
and the transverse impact parameters bi and b 2 , d'^(x/, x\) are the usual partonic cross sections for 
single quarkonium production and 5q^q^ is the Kronecker delta function. A further factorisation 
assumption is to decompose Yij{xi,X 2 , bi, b 2 ) into a longitudinal part and a transverse part 

Yij{xi,X2, bi, b 2 ) = Dij{xi,X2)Tij{\)u\)2), (2) 

where Dij{xi,X 2 ) is the double-parton distribution functions (dPDF) [47]. Moreover, by ignoring 
the correlations between partons produced from each hadrons, one can further assume 

Dij{XuX2) = fi{Xi)fjiX2), 

r,-/bi,b2) = r,(bi)rxb2), (3) 

where fi(xi) and fj(x 2 ) are the normal single PDFs. This yields to 

o-Q,Q. = -—\ -y. (Tik-^Q^o-ji^Q^ f d^b f T,{bi)Tk(bi -b)d^bi f Tj(b2)T/(b2 - b)d%. (4) 
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If one also ignores the parton flavour dependenee in Ti j^ iQi) and defines the overlapping funetion 


F(b) = J T(b^)T(b-b)d%, 

one reaehes the so-ealled “poeket formula” 


(5) 


(^QiQ2 


1 0'Qi0'Q2 

1 + ^(3i( 32 ^eff 


( 6 ) 


where o-q, and are the eross seetions for respeetively single (3i and Q 2 produetion and cTeff is 
a parameter to eharaeterise an effeetive spatial area of the parton-parton interaetions via 


O'eff - 



d^bF(bf 



(7) 


Under these assumptions, it is only related to the initial state and should be independent of the 
final state. However, the validation of its universality (proeess independenee as well as energy 
independence) and the factorisation in Eq.(6) should be cross checked case by case. In a fact, 
some factorisation-breaking effects have recently been identified (see e.g. [48, 49, 50]). Thanks 
to its larger luminosity and its probably wide rapidity coverage, AFTER®LHC provides a unique 
opportunity to probe DPS and to extract cTetf from double-quarkonium final states. 

To perform our predictions, we will use cr^ff = 5.0 + 2.75 mb, which was determined from 
7/i/f-pair production data at the Tevatron by DO collaboration [7].^ The reason for such a choice is 
that all of the double-quarkonium-production processes share the same gluon-gluon initial states 
and the typical x are not that much different. This also means that we only need to assume the 
energy independent of cTgff. However, we do not claim that this value is the only one possible; we 
only take it as our reference number. If one wants to use another value of CTgff, one can just simply 
perform a rescaling (proportional to l/cr^ff) of the numbers given in the following. 

Since the description of single heavy-quarkonium production at hadron colliders in the whole 
kinematical region is still a challenge to theorists, using ab initio theoretical computation of <tq 
would significantly inflate theoretical uncertainties. Instead, we will work in a data-driven way to 
determine ctq. 

Our procedure is as follows. We start from the cross section ctq. which can be written as 


o-{pp Q + X) 



dxi dX2fa(Xl )fb{X2) — \^ab^Q+X?dlAPSQ+x, 
2 s 


( 8 ) 


where fa, fb are the parton distribution functions (PDF) of the initial partons a and b, JLIPSq+x is 
the Lorentz-invariant phase-space measure for pp Q + X and is the partonic centre-of-mass 
energy (i.e. s = X1X2S). For single quarkonium production in p + p collisions at -^ = 115 GeV, 


^Note that Ref. [7] has updated the value of cTeff to be 4.8 + 2.55 mb. However, since the difference is very small, 
we still used the original one. 
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K 

A 

# of data 


J/lj/ 

0.67 + 0.08 

0.38 

51 

422 

(A(25) 

0.15 + 0.03 

0.35 

4 

1.12 


(a) Charmonia 



K 

A 

# of data 


Y(15) 

0.89 

0.084 + 0.0061 

14 

29 

Y(25) 

0.79 

0.056 

9 

2.2 

Y(35) 

0.68 ± 0.029 

0.046 

9 

3.9 


(b) Bottomonia 

Table 2; Results of a fit of cP'crIdPTdy to (a) the \p{nS) PHENIX data [51] by fixing n - 2 and {Pj) - 4.5 GeV and 
(b) the T(n5) data CDF [52] data by fixing n -2 and {Pj) - 13.5 GeV. Only the > 1% errors are given. 


the gluon-gluon initial state is dominant. The initial eolour and helieity averaged amplitude square 
for gg ^ Q + X ean be expressed in the form of a erystal ball funetion [20] 




K exp(-/c J) when Pt < (Pj) 

Q 

if exp(-/c^) (l + " when Pj > (Pt) 


(9) 


where K = A^ksIM^. The parameters K,A,n and (Pj) can be determined by fitting the (differen¬ 
tial) cross sections to the experimental data. The dedicated codes to perform the fit and to com¬ 
pute the DPS contributions to double-quarkonium production have been implemented in HELAC- 
Onia [53, 54]. 

Once a fit is done, is fixed and it allows us to evaluate cr(pp Q + X) (or its 

differential counterparts in any variable) which can then be injected into the “pocket formula” 
Eq. (6) in order to predict the DPS yield. Since we do not apply any muon cuts, we do not need to 
make any assumptions regarding the polarisation of the production quarkonia. 

The code was tested and, with the same parameters as in Ref. [20], we have reproduced their 
results. However, their combined fit of the charmonium data taken at the Tevatron and the LHC 
cannot reproduce well the low-energy data measured by PHENIX collaboration [51] at RHIC. 
Since the collision energy of RHIC ^/s = 200 GeV is very close to the centre-of-mass energy of 
the fixed-target experiment at the LHC (ALTER®LHC), i.e. = 115 GeV, we prefer to use the 
PHENIX data alone to determine the parameters in Eq. (9). A fit of (fcrldPTdy to the PHENIX 
data [51] for Jjifj and i/^(25) production gives the^f^ results presented in Tab. 2a having fixed n = 2 
and {Pt) = 4.5 GeV. We also show the comparisons of the Pj spectra in Eig. la-c. The large 
for the single 7/i/^ production can be reduced to 55.8 when one only considers the 23 PHENIX 
data points in the central region (i.e. < 0.35) and excluding the lowest-PT- bin. A fit to the 

sole PHENIX data in the forward region 1.2 < \yji^\ < 2.4 changes /c by ~ 15% and d by ~ 5%. 
However, the main uncertainty in predicting DPS contributions to double ij/ production remains 
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J/\j/ Comparison 


J/\|r Comparison 



\|/(2S) Comparison Y{1S) Comparison 



Y(2S) Comparison Y{3S) Comparison 



Figure 1: Comparisons with the PHENIX measurements [51] for 7/i/f (a,b) and t^(2S) (c) production and with the 
CDF measurements [52] for T(15) (d), T(25) (e) and T(35) (f) production. 

from that of and those from these fits are in praetiee nearly irrelevant for our predictions. This 
is obvious for A which only affects the normalisation. 

In contrast, there is no differential measurement of Y yields at RHIC. There exists data from 
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the fixed-target Fermilab experiment E866 [55] but only at low Pj- We therefore performed a fit 
of d^cr/dPrcly to CDF [52] Run I data at = 1.8 TeV. The results for Y are presented in Tab. 2b 
having fixed n = 2 and (Pt) = 13.5 GeV. For illustration, the comparisons between the fit and 
the CDF data [52] are shown in Fig. Id-f. Some comments about the fit are however in order. If 
we instead performed a combined fit to CDF [52], ATLAS [56], CMS [57] and LHCb [58, 59] 
data, the value of k (T) would be shifted by at most 30% (10%) but with significantly worse 
All this may however not be so relevant since, as for the charmonia, the fit to TeV data tend to 
underestimate the RHIC Rr-integrated Y production cross section as measured by STAR [60] by a 
factor a bit smaller than 2 - the STAR result has however a 30% uncertainty. The uncertainties on 
K and X given by the fit are therefore far too optimistic since the Crystall Ball parametrisation 
seems not to correctly capture the energy dependence of the cross section. The corresponding DPS 
yields of Y at AFTER®LHC which we give here should therefore be considered as conservative 
lower estimates. All of the above fits are performed with MSTW2008NLO PDE set [61] avail¬ 
able in LHAPDP5 [62] and the factorisation scale Hf = -l- R|. The physical mass Mq for 

quarkonium is taken from PDG data [63] as well as the branching ratios. 

2.2. Single-parton scatterings 

2.2.1. Double-charmonium and double-bottomonium production 

The SPS contribution to 7/(/f-pair production have systematically been investigated in our pre¬ 
vious works [16, 18]. We have shown that a leading order (LO) calculation in the strong coupling 
constant, as, is enough to account for the low-Rr data as well as the Rr-integrated cross section, 
the bulk of the events lying at low Pt- However, if one goes to mid Pt (e.g. Rj- > 5 GeV), 
0(al) contribution start to be large. As a consequence, the yield and the polarisation changes 
significantly compared to a LO calculation. Since we are only interested in the data which are 
measurable with up to 20 fb“* in order to assess the feasibility of measuring quarkonium-pair pro¬ 
duction with ALTER@LHC, we will focus on the low Pt region. As we will explicitly show, LO 
evaluations happen to be sufficient. Besides, the colour-octet contributions are also negligible at 
low Pt for they are suppressed by powers of v without any kinematical enhancement at variance 
with the single-quarkonium-production case. 


decay channel branching ratio (%) 


<A(25)- 

Jlilf + X 

57.4 

Y(2S) - 

> Y(1S)+X 

30.2 

Y(35) - 

> Y(1S)+X 

8.92 

Y(3S) - 

> Y(2S) + A 

10.6 


(a) Decay within 

a family 


decay channel 

branching ratio (%) 

7/(A ^ p^p~ 

5.93 

1^(25) p^p~ 

0.75 

rilS)^p^p- 

2.48 

ri2S)^p^p- 

1.93 

Y(35) ^ juV 

2.18 


(b) Leptonic decays 


Table 3: Various decays (and branching ratios) considered in this article [63]. 

On the contrary, the feed-down contributions from higher excited spin-triplet S -wave quarko¬ 
nium has to be considered. It is substantial as already shown for the Jli/f-pair production in 
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Ref. [18]. These will systematieally be taken into aeeount in our predietions as done in Ref. [18]. 
The branehing ratios that will be used in this eontext are taken from PDG [63] and we have listed 
them in Tab. 3 for completeness. 

The general formula for the amplitude of the production of a pair of colour-singlet (CS) S -wave 
quarkonia <3i and Q 2 with as initial partons a and b is 


•^ab^Q'l\Pi)+Q2^{P2)+X 


( 10 ) 


E E 

il,i2,Cl,C2 i3.i4.C3.C4 


A^(/li|.yi,.y2)A^(d2|.y3.'y4) dc|C2<^C3C4 ^l(0)^2(0) 

^Mq,Mq, N, An 


e)^(Pi=o)+Q)]eJ](P2=o)+x’ 


where we denote the momenta of quarkonia (3i and <32 as Pi and P 2 respectively and their polari¬ 
sations as Ti 2 , N(Ai^ 2 \si, 3 , ■^ 2 , 4 ) are the two spin projectors and Pi, 2 ( 0 ) are the radial wave functions 
at the origin in the configuration space for both quarkonia. In the above equation, we have defined 
the heavy-quark momenta to be < 2 'i, 2 , 3,4 such that Pi 2 = ^ 1,3 + ^ 2,4 and pi ^2 = (?i ,3 - ^ 2 , 4 )/ 2 . i'i, 2 , 3,4 
are then the heavy-quark spin components and dacj/ is the colour projector. The spin-triplet 
projector N(A\si, sj) has, in the non-relativistic limit, v ^ 0, the following expression 

N{A\si, sj) = —p— Si). (11) 

2 V2Mq 2 2 

All these computations can be performed automatically in the HELAC-Onia [53] framework based 
on recursion relations. The radial wave functions at the origin P(0) are taken from Ref. [64], which 
were derived in the QCD-motivated Buchmilller-Tye potential [65]. We also listed their values in 
Tab. 4. 


Quarkonium 

|R(0)|2 (CeVQ 

■//<A 

0.81 

«A(25) 

0.529 

T(15) 

6.477 

T(25) 

3.234 

T(35) 

2.474 


Table 4; The radial wave functions at the origin squared |/?(0)|^ [64] of S -wave quarkonium used in this article. 


2.2.2. Charmonium-bottomonium pair production 

The simultaneous production of a charmonium and a bottomonium has been studied in Refs. [13, 
19]. Its CSM contributions are expected to be suppressed because the direct LO contributions in 
CS mechanism (CSM) are i.e. suppressed compared to double-charmonium and double- 
bottomonium production. Hence, it is expected to be a golden channel to probe colour-octet 
mechanism (COM) at the LHC [13]. However, such a statement is valid only if one can clearly 
separate DPS and SPS events experimentally since the DPS contributions would be substantial. 
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For a thorough discussion, the reader is guided to [19]. In eontrast, eolour oetet (CO) eontribu- 
tions ean appear at whieh however are suppressed by the small size of the CO long distanee 

matrix elements (LDMEs). If one follows the arguments of Ref. [13], one is entitled to consider 
only the cc(^ + bb^S ®), ccCs + bb{^ and cc(^ + bb(^ ehannels. This approxima¬ 
tion is however based on the validity of the veloeity sealing rules of the LMDEs whieh may not be 
reliable. A complete eomputation -even at LHC energies- aecounting for all the possible ehannels 
up to in NRQCD is still lacking in the literature: there are indeed more than 50 ehannels at LO 
in as eontributing io ip + T produetion. Thanks to the the automation of HELAC-Onia [53, 54], 
sueh a eomplete is at reaeh. 

The formula for the S -wave CO amplitude is similar to that for CS state produetion with the 
following formal replaeements for CO in Eq.(ll) 




CiCj 


Ri(0) 

V4^ 


V(27+1)(A2-1)’ 


( 12 ) 


where T^.^. is the Gell-Mann matrix and (C)'(^®)) is the CO LOME. We refer the reader to 
Ref. [53] for the E-wave amplitudes. 

The non-perturbative CO LDMEs should be determined from experimental data. Their values 
unfortunately depend mueh on the fit proeedures. We took four sets of LDMEs from the literature 
(see the details in Appendix A. 2). 

Einally, we deseribe our parameters for our SPS ealeulations. In the non-relativistie limit, the 
mass of the heavy quarkonium ean be expressed as the sum of the eorresponding heavy-quark-pair 
masses. In our ease, we have 

Mq = Imq, (13) 

where mg = for eharmonium and mq = for bottomonium. The masses of eharm quark 
and bottom quark are taken as = 1.5 + 0.1 GeV and nib = 4.75 + 0.25 GeV. The fac¬ 
torisation scale jjiF and the renormalisation scale ijlr are taken as fip = bR ^ with 

yUo = The advantage of using /io = is that we are 

able to reeover the eorreet mass threshold MQ^ + Mq^ in the low Pj regime. Einally, the PDE set 
for the SPS ealeulation is CTEQ6L1 [66] with the one-loop renormalisation group running of as- 


3. Energy dependence of the ratio DPS over SPS 

Due to the very large integrated luminosity of ALTER®LHC (up to 20 fb“' per year) eom- 
pared to the experiments performed at RHIC, the measurement of double-quarkonium produetion 
at AETER@LHC will provide a unique test of the interplay between the DPS and SPS produc¬ 
tion mechanisms in a new energy range. The energy dependenee of CTeff will be explored at a 
wide energy range when eombined with the LHC collider and Tevatron data^. Due to the double 


^ Since we noted that the energy dependence obtained with the partonic amplitude {gg —> QX) given by a Crystal 
Ball fit with fixed parameters is not optimal when going to TeV energies down to RHIC energies, we have used the fit 
parameters of [20] (based on a fit of Tevatron and LHC data) to predict the DPS yield in the TeV range and our fit to 
the PHENIX data for the RHIC and fixed-target-experiment energy range. 
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Vs [GeV] 


Figure 2: (Upper panel) The cross sections of (prompt-)7/i/f pair production via SPS and DPS mechanisms for two 
values of cTeff as a function of (Lower panel) DPS over SPS yield ratio for 5 < CTeff <15 mb. The black circles 
correspond to 10 mb. [Aside from the choice of cTsS, no theoretical uncertainties are included]. 


enhancement of the initial gluon-gluon luminosity with the energy, ^fs, DPS contributions are 
expected to be more and more important with respect to the SPS ones at larger ^fs. This can be 
observed on Fig. 2. 

One however sees on Fig. 2 that a change of (Teff from 15 mb -which seems to be the favoured 
value for jet-related observables- to 5mb -which is the value extracted by DO from the 7/i/f -l- J/ij/ 
data [7]- results in a significant change in the point where both contributions are equal. In the 
former case, it occurs very close to the energy of AFTER@LHC, in the latter case, it occurs 
between the Tevatron and the LHC energies. All this clearly motivates for measurement and cre^ 
extractions at low energies. 

4. Impact of the QCD corrections at low transverse momenta 

Before showing our results and in order to motivate the use of LO predictions for this ex¬ 
ploratory study, we have found it useful to give an explicit comparison between the differential 
cross section at LO and NLO* for double-7/i/f production in the kinematical domain accessible 
with 20 fb“\ that is up to transverse momenta on the order of 10 GeV at the very most. Indeed, 
in a previous study [16], we have showed that the impact of the real-emission corrections, such as 
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I Ay I between ijj 



(a) Absolute rapidity difference between both 7/^/r 


(b) Pair Invariant mass 



Leading Pj of two ip (GeV) 



(c) Leading Pj among the 7/i/r pair 


(d) Pair transverse momentum 


Figure 3: LO vs. NLO* differential distributions. 

gg JI iff + J/i^ + g, becomes increasingly important at large transverse momenta. 

Figs. 3 show the comparison between LO results and NLO* results (which are known to repro¬ 
duce well the full NLO [17]). The invariant-mass and rapidity-difference spectra are not affected 
by the real emission at . Indeed, in the low-Fr region, the Born topologies are dominant, and 
there is no kinematical enhancement in the real-emission topologies which could compensate the 
as suppression. Only when one goes to large transverse momenta, these are enhanced and can 
become dominant. This explains the difference in the slope as a function of the leading Pt in 
Fig. 3c. The results are however similar for Pj < 10 GeV where the cross sections are larger than 
0.1 fb. 

In addition, as we already discussed in Ref. [18], at LO, a 2 ^ 2 ki nematics for SPS would 
result in a transverse momentum of the 7/i/^-pair P'lf being zero and in a trivial LO distribution 
on Fig. 3d. This is however not the case if one takes into account a possible intrinsic kj of the 
initial partons which can also been considered as a part of QCD radiative corrections - initial- 
state radiations to be precise. Such a smearing can be phenomenologically be accounted for and 
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compared to a pQCD result. To do so, we have smeared the kinematics of LO events using a 
Gaussian distribution with {kj) = 1 & 2 GeV as done in Refs. [16, 18]. We stress that the value of 
{kj) is essentially empirical, hence the choice of two values for illustration (resp. curves labelled 
sml and sm2). This can thus be compared with our NLO* curves in the accessible domain with 
C)(20)fb“^ at AFTER@LHC, that is < 10 GeV. One sees that the smearing mimics relatively 
well the effect of the QCD corrections with {kj) = 2 GeV which we will use in the following for 
the comparison with the DPS yield. Overall, the P^j distribution is obviously very different than 
a single peak at 0. 

5. Predictions at AFTER@LHC 

We are now in the position to present our numerical results at ^fs = 115 GeV in p+p collisions. 
The total cross section we obtained are given in Tab. 5, 6 and 7. The results have been multiplied 
by the branching ratios into a muon pair and they are all in unit of fb. In general, we have 

The DPS contributions decrease quickly when the mass threshold Mq^ + Mq^ increases because 
of its square dependence of the initial-state parton luminosity. With the nominal integrated lu¬ 
minosity of 20 fb“' proposed to be collected at AFTER@LHC, we find that the measurement 
double-bottomonium production is out of reach‘d and one may be able to record a few J /ip + 'T{IS) 
events, which receives substantial DPS contributions. 



7/0- + J/ip 

J/ip + ip(2S) 

ipi2S) + ipi2S) 

O' DPS 

^CSM 
^ SPS 

590 :™ 

700!“» 

19+23 

oc+440 

°^-68 

^■'--’-0.052 

2 ^+13 


Table 5: cr(pp Qi + Q2 + X) x !B(Qi — » p^p )S(d ?2 ^ p^p ) in units of fb at 2 fs - 115 GeV, where <3i,<32 = 
Jlp, p{2S). The DPS uncertainties are from cr^ff and the SPS ones from niQ and the scales. 


One should however always keep in mind that ctsps for i/r -i- T production strongly depends 
on the CO EDMEs. We have investigated this dependence in Appendix A.2 with four different 
sets of EDMEs and the results vary up to one order of magnitude which precludes any strong con¬ 
clusions^. In addition, these EDMEs are usually fit from the experimental data at high transverse 
momentum region and are known to overestimate the single-quarkonium yields at low Pj (see [67] 
and references therein). This is also probably the case for quarkonium-pair production especially 
when they come from single gluon splittings. We have therefore find it only meaningful to show 
upper limits on ctsps for ^ -b T production in Table. 6. These numbers are in any case at the limit 
of observability. 


^We note that such a measurement has never been done anywhere else. 

^For convenience and possible future studies, we have tabulated in Appendix A. 1 the values of all the relevant 
short-distance coefficients which can then be combined with any LDME set. 
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7/(A + T(lS) 

7/(/( + T(25) 

7/(A + T(3S) 

O" DPS 
NRQCD 
•^SPS 

0 17+0-21 

-0.058 

<0.69 

0 037'^° ‘)45 

D-DO /_Q QJ3 

< 0.14 

0-018t°^3^ 

< 0.11 


(/r(25) + T(15) 

(A(25) +Y(25) 

(/((25) + T(35) 

O" DPS 
NRQCD 
•^SPS 

2 6 • 10-3 +3.2-10-3 

^ -9.1-10-4 

< 0.031 

5 7 10“^ +6.9-10-4 

-2.0-10-4 

< 5.4 ■ 10-3 

7 S 1 n“4 +3.4-10-4 

2-8 • 10 _9g.jQ-5 

< 3.0 ■ 10-3 


Table 6: cr(pp —> (3i + (32 + ^) x S((3i ^ p^p~)S(Q .2 —> in units of fb with -^ = 115 GeV, where 

(3i = Jlif/,tf/{2S) and Qi - T(15), T(25), T(35). For SPS production, only the upper limits of the yields are given 
(see text). The DPS uncertainties are from a-^g. 


The quoted theoretieal uneertainties in the tables result from the variation of (Tgff within 5+2.75 
mb for the DPS yields and from the scale uncertainties as well as heavy-quark-mass uncertainties 
for the SPS yields, as discussed in Sec. 2. 

As regards double-charmonium production, about 10 thousand events could be collected per 
year -which is more than what has so far been collected by LHCb and CMS. In the analysis 
of the differential distributions, we therefore only focus on these and, in particular, on 7/i/r-pair 
production. We show three interesting distributions without kinematical cuts. Along the lines 
of [40], we also used the LHCb kinematical acceptance, i.e. the rapidity of Jjif/ restricted to be in 
the interval of [2,5]. 



T(15) +Y(15) 

T(25) + T(25) 

Y(35) +Y(35) 

O' DPS 

1-2 • 

5 6 ■ 10-’ +6.8-10-3 

1 4.10-2 +1-7-10-' 
-4.7-10-8 

^CSM 
^ SPS 

2.8 • 10-3 :^-3;;o- 

3 5 1 n“4l +1.7-10-3 

-2.8-10-4 

9 9 1 n“4l +1.1-10-3 

-1.8-10-4 


T(15) + T(25) 

Y(15) + T(35) 

Y(25) +Y(35) 

O' DPS 

c -1 -| 0“® +6.2-10-® 

^ -1.7-10-6 

9 5 1 n“6 +3.0-10-6 
-8.7-10-3 

5 5- 10-2 +6-7-10-3 

^CSM 
^ SPS 

2 0 • 10-3 +9-3-10-^ 

^ -1.6-10-3 

1 6 ■ 10-3 +7.4-10-3 

^ -1.3-10-3 

5.6 • 10-^ +2.6.10-3 

-4.4-10-4 


Table 7: cripp — » (3i -H (32 + X) x S((3i — » p^p )S(Q2 — » p^p ) in units of fb with -\/7 = 115 GeV, where 
(3 i,( 32 = T(15), T(25), T(35). The DPS uncertainties are from cr^g and the SPS ones from the mg and the scales. 


The absolute rapidity difference between the JjiJ/ pair is expected to be a good observable to 
discriminate the DPS and SPS contributions. This was first pointed out in Ref. [20] and this was 
used later on by DO collaboration [7] to extract crgg from double-7/(/r production at the Tevatron. 
The DPS events should have a broader distribution in Ay than the SPS ones, because two (rela¬ 
tively) independent hard interactions happen simultaneously in DPS while the two J/if/ from SPS 
are more correlated. The situation still does not change at AFTER@LHC without or with cut as 
Fig. 4 (left) and (right) show. In the latter case, the restriction to negative rapidities in the centre- 
of-mass obviously reduce the Ay range. Starting from Ay = 2, the DPS events dominate the SPS 
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Figure 4: Differential cross section as a function of the absolute rapidity difference of the Jlt/j pair, without (left) or 
with (right) a rapidity cut. 




Figure 5: Differential cross section as a function of the invariant mass of the Jjtj/ pair, without (left) or with (right) a 
rapidity cut. 

events. A ratio DPS/SPS of 10 is obtained for Ay > 2. The distribution of the invariant mass for 
the J/ij/ pair refleets a similar information as the Ay distribution. Hence, it follows that the 
spectra of DPS are also broader than those of SPS, which can be seen on Fig. 5 (left) and 
(right). 

As we discussed earlier, predictions for the dependence of the SPS yield depend much 
on the kr smearing of the initial partons which can mimic a part of the QCD corrections. Due to 
the relative smaller yields at AFTER@LHC energies than at LHC energies, one can only access 
PP' < 10 GeV, as illustrated on Fig. 6. In such a kinematical region, the kr smearing effect makes 
the SPS spectrum as broad as the DPS one with (kr) = 2 GeV. 

Finally, we present on Fig. 7 the cross section as a function of the total rapidity of the J/if/ pair 
(right), and of the sub-leading Pj between the Jjif/ pair (left). One sees that the sub-leading Pj 
spectrum may be measured up to 6 GeV with AFTER@LHC. As regards the rapidity distribution, 
its maximum is obviously located at Fcms = that is F = 4.8 in the laboratory frame. One sees that 
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Figure 6: Differential cross section as a function of the transverse momentum of the Jjtj/ pair, without (left) or with 
(right) a rapidity cut. 



Figure 7; Differential cross section as a function of (left) the sub-leading Pj with a rapidity cut and (right)the rapidity 
of the 7/i/f pair. 

one can expect some counts down to ^ 2.5 where xp - sinh(y^^ - 4.8) - -0.5. This is 
precisely the kinematical region where double intrinsic cc coalescence contributes on average [10]. 
Any modulation in the pair-rapidity distribution would sign the presence of such a contribution. 

Finally, we have investigated the impact of using different (double)PDFs (MSTW2008NLO [61], 
CTEQ6L1 [66], GS09 dPDF [47]) on differential distributions are also shown in Fig. 8; they are 
found to be moderate in all cases. 

6. Conclusion 

We have discussed double-quarkonium production in proton-proton collisions at a fixed-target 
experiment using the LHC proton beams, AFTER@LHC. These processes have lately attracted 
much attention, both in the theorist and experimentalist communities. They are expected to be 
good observables to further constrain the various models describing heavy-quarkonium produc¬ 
tion. Double-quarkonium production also provides a good opportunity to study DPS since the 
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(a) Pair transverse momentum 


(b) Absolute rapidity difference between both Jl^j/ 
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Figure 8: Differential distributions for DPS with various PDFs: (a) transverse momentum spectrum; (b) absolute 
rapidity difference ; (c) invariant mass distribution; (d) rapidity of Jj^ pair. 

yields of single quarkonium produetion is large and their deeay to four muons is a elean signal at a 
hadron eolliders. AFTER®LHC provides very appealing opportunities to study these observables 
with a LHCb-like deteetor and in new energy region. 

In this paper, we have studied both DPS and SPS eontributions for double-quarkonium pro- 
duetion. These proeesses inelude ipiriiS) + if/(n 2 S), iJ/{niS) + T(mi5) and Y(mi5) + T(m25) with 
ni,n 2 = 1,2 and mi, m 2 = 1,2,3. DPS contributions are estimated in a data-driven way, while SPS 
ones are calculated at LO in non-relativistic QCD (NRQCD) [68], more precisely in the CSM for 
il/{niS)-\-ilj{n 2 S) and Y(mi5)-l-Y(m25) and accounting for CO contributions for -l-Y(mi5). 

From our calculations, we find that ten thousand of double-charmonium events can indeed be mea¬ 
sured at AFTER®LHC with the yearly integrated luminosity of 20 fb“k In the most backward 
region, a careful analysis of the rapidity distribution could also uncover double intrinsic cc coa¬ 
lescence contributions. In general, future measurements on double-charmonium production can 
provide extremely valuable information on QCD, in particular important tests on the factorisation 
formula for DPS and the energy (in)dependence of crgif. 
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Appendix A. Charmonium-bottomonium pair production in NRQCD 

Appendix A.l. Short-distance coefficients for charmonium-bottomonium pair production 

In NRQCD [68], the cross section for a charmonium C and a bottomonium S production can 
systematically be written as 

cr(C + S) = ^^(T(cc[ni]bb[n 2 ]) X {Cfffii)) X {0^{n2)), (A.l) 

n\,n2 


where ni,n 2 are different possible Fock states, cr(cc[ni] + bb\n 2 \) is the short-distance coefficient 
(SDC) for the production of a charm-quark pair in the Fock state ni and a bottom-quark pair in 
the Fock state n 2 simultaneously. The LDMEs {Cfffii)) and {Cffi{n 2 )) should obey the velocity¬ 
scaling rules of NRQCD. In this appendix, we give the complete list of the SDCs for charmonium- 
bottomonium pair production at the leading order in in proton-proton collisions at the center- 
of-mass energy of yfs = 115 GeV. This includes the contributions from = 

0,1,2) for S -wave quarkonium production and from ^ = 0,1,2) for F-wave quarkonium 

production. There are in total 66 non-vanishing channels to be computed. Such a computation is 
automatic in HELAC-Onia [53, 54], but has never been carried out even at LHC energies. Thanks 
to the heavy-quark-spin symmetry of NRQCD, we have 


^^,S(3p[8])) ^ ^27 + 1) X 


(A.2) 


We can thus define new SDCs relevant for P' 


3n[8] 


cr(cc[^ F™] -I- bb[n 2 ]) = ^ (27 -l- 1) x cr(cc[R™] -l- bb[n 2 ]), 

7=0 7=0 

2 2 

cr(cc[ni] -I- bb[^ ^R®]) = ^ (27 -i- 1) x cr(cc[ni] -l- hh[^®]). 


(A.3) 


7=0 


7=0 


Therefore, we have 

2 


^ m(cc['R®] + bb[n2]) X <0^('R®))<0V2)) = fr(cc[^ + bb[n2]) X <0^('R™))<0®(n2)), 


7=0 

2 


7=0 


^ m(cc[ni] + bb[¥f^]) X {OP(n,)){0^(¥f^)) = (r(cc[n,] + bb[J] x {OP(n,)){<¥(¥^^)). 


7=0 


7=0 


(A.4) 


We display the numerical values for the SDCs for these Fock states in Table. A.8 with CTEQ6L1 [66] 
as our PDE set. 
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Fock state 

bb[¥l^] 


bb[¥^^] 

biH.PXT'i 

bb[¥^^^] 

bb[¥l^] 

bb[¥^^^] 


- 

i3;« 

- 

- 

- 

- 

- 

CC[^®] 

cc[!s®] 

cclH”] 

cc[>f/^] 

40+200 

770+4000 
' ' -620 
770+1100 
-170 

470+2700 

■31+180 
-’ + -25 
91+130 
^ + -18 
91+120 

4' + -18 

7700+14000 

z, / uu_22oo 

i2oo;g“ 

2io;ir 

69;«« 

7 r+310 
/-J-40 

770+4200 
' -590 

180!”«o 

oon+2400 

'^'^^-280 

7 C+I 8 O 

^^-21 

7 5 +^"^ 

-^-6.4 

+’■+-5.2 

160«» 

46:g" 

01 +220 
-’ + -26.0 
12+87 
+^-10 

3.6:r, 

3 . 0 X 

7 0+44.0 
^•'^-6.0 

1.2+f6 

0 37+2-^ 

-0.31 

0 33+^-^ 
u.-1-1-028 

0 15+^-^ 

-^-0.13 

43,go 

Q 1 +56 
-7.6 

Q +58 
°'-6.8 

-1.1-2 7 

+ • -0.86 

0 79+^-^ 
+’■'’^-0.68 


Table A.8; The SDCs (at the leading order in a^) for the various combinations of the Fock states contribut¬ 
ing to charmonium-bottomonium pair production at -v/i = 115 GeV. The unit of the SDCs of cc[ni] + bb[n2] is 
fb/Gewhere L, = 0 when n, is 5-wave and L, = 1 when n, is F-wave. The uncertainty quoted is coming 

from the variation of jif - jir e [5/io,2yL(o] (pLo — + irih)^ + P\) and the uncertainties on iric - 1.5 + 0.1 GeV 

and nib = 4.75 + 0.25 GeV. 


Appendix A.2. Cross sections for single-parton scattering 

From the SDCs given in Table. A.8 and the LDMEs extraeted from the experimental data, 
we are now able to estimate the eross seetions of eharmonium-i-bottomonium pair produetion at 
= 115 GeV. The values of the LDMEs however signifieantly differ depending on the different 
experimental input data and the different fit setup. Eor example, the CO LDMEs of J/f extraeted 
from pp data can be quite different with or without NLO QCD corrections. Here, we will discuss 
the results based on four sets of LDMEs for charmonia and bottomonia, which can be described 
as follows: 

Set I: This set is based on the LDMEs of J/tff,if/(2S) dxxd Xc presented in Ref. [69] and those of 
T(1S), T(2S), T(35) S-nd Xb{^P),Xb{'^P) presented in Ref. [70]. They are extracted from 
Tevatron data with SDCs at LO in Us- The LDMEs of XbOP) have been set to zero in the fit 
of Ref. [70] ^ 

Set II: This set is based on LDMEs of 7/(A,iA(25),Trc,T(I5),T(25),T(35),Trfa(lR),T^>(2E) pre¬ 
sented in Ref. [72]. The contributions of XbOP) have been ignored. Hence, we will set 
the LDMEs of XbOP) to be zero. The fit was performed at LO in The LHC, Tevatron 
and RHIC data were used to perform this combined fit. 

Set III: This set is based on LDMEs extracted from NLO analyses, i.e. the LDMEs of J/f, f{2S ),Xc 
from Ref. [73] and those of T(nS),Xb(nP),n = 1,2,3 from Ref. [74]. The CO LDMEs of 
charmonium are extracted from Tevatron data [73], while both Tevatron data and LHC data 
were used in Ref. [74]. 

Set IV: This set is based on LDMEs for charmonium [75] and bottomonium [76] production based 
on other NLO analyses. They are determined by a combined fit to Tevatron and LHC data. 


®Note that both fits used CTEQ5L [71] whereas we have used here CTEQ6L1, whose results are anyhow very 
close. 
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j/t^ + r(is) 

J/i// + r(2S) 

j/i// + r(3S) 

Set I 


0.0185!S;™2 

0.0158:S 

Set II 

0 0948+°-^^' 
u.u:?‘to_Q Qg26 

0.0146!SS 

6.28 ■ 10-3 

Set III 

u.u / u / _Q Qg75 

0.0205!S;;15, 

1,14.10-2+6.34.10-2 

Set IV 

-0.0163 

6 00 • 10“^ +3.3610-2 

U.UU iU _4g9,io-3 

2 51 • 10-3 +1.3410-2 
Z,.Ji iU _2 03.10-3 


iA(2S) + T(15) 

|A(2S) + T(2S) 

iA(25) +Y(35) 

Set I 

1 8S 1 +1.0110“^ 

1.85 lU _i 50.10-3 

5.83 • 10-^ 

4 64 • 10-4 +2.57.10-2 

‘+.U‘+ ru _37g.io-4 

Set II 

4 30 • 10“^ +2.6210-2 
H-.JV iU _373 .io-3 

6-78 • 10-^ 

7 OQ • 10“'^ +1.64.10-2 
j.vy iu _2 49.10-4 

Set III 

3 19- in~2 +1.9810-2 

J.Ly iU _2g4.io-3 

817.10-4+4.^.10-3 

4 57 • 10-4 +2.54.10-2 
-4.1110-4 

Set IV 

9 03 • 10“^ +4.78-10-2 

y.vj iu _730 .io-4 

2.80 • 10-4 

1 49 . lO""^ +6.81.10-4 

iU _ii3.io-4 


Table A.9: crspsipp —> (3i + Q2) x !B(Qi —> p^p )S(<32 —> P^P ) in units of fb with a/s = 115 GeV, where 
Qi = Jli//,il/(2S) and Q 2 = T(15), Y(25), Y(35). We take four sets of LDMEs. 



7M +Y(15) 

j/i// + r(2S) 

7M +Y(35) 

{^[l]^3^[8]j 

6 1 • 10-^ + 3 . 010-2 

-4.910-2 

1 8 10-3 +8.610-2 

1.8-10 _i4.io-3 

2 5 - 10-^ + 1 . 2 . 10-2 

exclude feeddown 

0 024'^°45 

u.uz,‘t-0020 

6 0 • 10-^ +3.710-2 

^-^^^-8.810-2 

include feeddown 

0 060+°-^^ 

0 019‘^°41 

U.ui5'-0015 

0 016'^°‘)95 

o.oiu_QQj3 


iA(2S) + T(lS) 

(A(25) +Y(2S) 

|A(2S) + T(35) 

{^[l]^^[8]j 

6 1 • 10-4 + 3 . 010 -^ 

-4.910-4 

1.8 • 10-4 +^o f , 

-l.ScdotlO ^ 

2 4.10-4 +1.2-10-2 

-1.910-4 

exclude feeddown 

1 1 • 10-3 +6.110-2 

• ^ -9.110-4 

3 0 • 10-4 +1.610-^ 

2’-0 40 _24.io-4 

4 6 - 10-4 +2.610-2 

^ -3.810-4 

include feeddown 

1 9 • 10“2 +1.010-2 
4-^ -1.510-2 

^8 10-4 +3.2.10-2 

5.8-10 _4 7.io-4 

4 6 - 10-4 +2.610-2 

-3.810-4 


Table A.10: cr^psipp —> (3i + Q 2 ) x !B(Qi —> p^p^)!B(Q 2 —> p^p^) in units of fb with V* = 115 GeV, where 
(3i = Jlif/,i^(2S) and Q 2 — Y(15), Y(25), Y(35). We have used the Set I of the LDMEs. The uncertainty quoted 
comes only from the SDCs. 
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In order to take into aeeount the feeddown eontributions, we have taken the neeessary braneh- 
ing ratios from PDG [63] . For the unknown branehing ratios, sueh as 'QrixbOP) T{nS ) + y), we 
used the estimated values from Table I of Ref. [74]. The SPS eross sections of ^ + T production in 
proton-proton collisions at = 115 GeV are presented in Table. ??. As clearly shown, the cross 
sections significantly differ from one set of LDMEs to another. Before closing this appendix, we 
would like to stress several points. 

• Because some CO LDMEs in Set II and Set IV are negative, the cross sections might be 
negative, which is of course unphysical. Eor example, the cross section for direct J/ifj -l- 
T(25) production (which then excludes feeddowns) is negative for the Set II and Set IV. 

• If one follows the arguments of Ref. [13], one is entitled to consider only the cc(^®) -l- 

(?&[^®], cc[^j^^] -I- (?&[^®] and cc[^®] -i- (?5[^ j^^] channels. This approximation is how¬ 
ever based on the validity of the velocity-scaling rules of the LMDEs which may not be 
reliable. By using Set I of LDMEs, we have shown the comparison in Table. A. 10. The row 
^ ^ ® only include cc^S f^)+bb[^ f’], cc[^ [^]+bb[^ ®] and cc[^ f^]+bb[^ J^’] channels, 

while the remaining lines contain all CO and CS contributions (with or without feeddown 
contributions). The results clearly show that the cc^S + bbfS ®], cc[^ j^^] -l- bb[^ ™] and 
cc[^®] -I- (7^[^j^^] channels are not sufficient. Moreover, the feeddown contributions are 
also significant but for ij/{2S) + T(3S). 

• The CO LDMEs used in this section are mainly determined by data in the high transverse 
momentum region. It is important to point out that these LDMEs yield to cross sections 
overestimating the data in the low transverse momentum region and, hence, the total cross 
sections for the single quarkonium production (see e.g. Ref. [67]). Hence, it is likely that 
any such NRQCD based estimation of -l- T at low Pj are too optimistic. However, as a 
conservative estimation, it is reasonable that we consider them as conservation upper limits 
of the SPS contributions (see Table. 6). 

• Einally, let us note that the relative importance of pure CO-i-CO contributions as compared 
to the mixed CO-i-CS depends much on the LDME sets. It essentially ranges from 30 to 
70 % irrespective of the charmonium-bottomonium pair which is considered. Eor the sake 
of completeness, let us add that the pure CS-i-CS from double feed-down from;^^ +Xb is on 
the order of a couple of per cent, but for Set IV where it can be up to 10%. 
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